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A DOUBLET LATTICE METHOD FOR THE DETERMINATION 
OF ROTOR INDUCED EMPENNAGE VIBRATION AIRLOADS - 
ANALYSIS DESCRIPTION AND PROGRAM DOCUMENTATION* 


SUMMARY 


An efficient state-of-the-art method has been developed to determine 
the unsteady vibratory airloads produced by the interaction of the main 
rotor wake with a helicopter empennage. This method has been incorporated 
into a computer program. Rotor Induced Empennage Vibration Analysis (RIEVA) . 
The program requires the main rotor wake position and the strength of the 
vortices located near the empennage surfaces. A nonlinear lifting surface 
analysis is utilized to predict the aerodynamic loads on the empennage 
surfaces in the presence of these concentrated vortices. The analysis 
has been formulated to include all pertinent effects such as suction of the 
interacting vortices and the shed vorticity behind the empennage surfaces. 

The analysis employs a stepwise solution (time domain); that is, a period 
corresponding to one blade passage is divided into a large number of time 
intervals and unsteady airloads are computed at each step- The output of 
the program consists of chordwise and spanwise airload distributions on 
the empennage surfaces. The airload distributions are harmonically analyzed 
and formulated for input into the Coupled Rotor/Airframe Vibration Analysis 
(Ref. 6). This report describes the theoretical development for the analysis, 
the RIEVA program documentation and output of a sample case of the program. 


*The research effort which led to the results in this report was financially 
supported by the Structures Laboratory, USARTL, (AVRADCOM) . 


I • INTRODUCTION 


The unsteady airloads produced by the interaction of the main rotor wake 
with the helicopter empennage (horizontal stabilizer/vertical fin) can be a 
major source of vibratory loads on a helicopter (Ref. 1). An efficient state- 
of-the-art method has been developed to determine these unsteady airloads due 
to the passage of concentrated main rotor wake vortices in the vicinity of 
the empennage. Since future helicopters are expected to be operating with 
high disk loading resulting in strong blade tip vortices in the rotor wake, 
the order of magnitude of these unsteady airloads will be significantly 
higher. As a consequence, these airloads will have to be determined in the 
early design and development stage to insure an efficient design. Therefore, 
the main purpose behind the development of this analysis is to establish an 
effective reliable methodology which can be utilized to analytically predict 
these rotor wake induced empennage airloads. 

The prediction of these aerodynamic forces requires suitable analysis 
techniques for the description of the nonuniform flow environment in which 
the empennage operates. The technique used herein involves the utilization 
of two basic programs to study these main rotor wake/empennage interaction 
phenomena, A wake analysis (Ref. 7) is used to determine the main rotor 
wake position and the strength of the vortices that are located near the 
empennage surfaces, and a nonlinear lifting surface analysis is utilized to 
predict the aerodynamic loads on empennage surfaces in the presence of these 
concentrated vortices. This nonlinear analysis, called RIEVA (Rotor Induced 
Empennage Vibration Airloads), has been formulated to include all pertinent 
flow phenomena such as the suction effects of the interacting vortices, span- 
wise induced effects, etc- For the steady-state flow case, the basic method 
used has predicted accurately the pressures on a swept wing in the presence 
of prescribed vortices as indicated by the results shown in Fig. 1 (Ref. 2 
and 3). The present effort involves the extension of this analysis to account 
for the unsteady main rotor wake induced effects and also the effect of shed 
vorticity behind the empennage surfaces. The present method is an improve- 
ment over recently developed methods (e.g., Refs. 4 and 5) for the prediction 
of airloads because it includes nonlinear effects due to the interacting main 
rotor tip vortices. 
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LIST OF SYMBOLS 


a constant, 


1/27t^. 


an element of the influence matrix. 

mean chord length of lifting surface, feet . 

lift, pitching moment (about quarter-chord) and rolling 
moment coefficients, respectively, of the lifting surface 
nondimensionalized to (q * Area), (q * Area * c) and 
(q * Area * c) respectively, where q = l/2pV^ , 

lift, pitching moment and rolling moment coefficients, 
respectively, corresponding to total loads (potential 
plus suction airloads) . 

2 

pressure coefficient, 2 (po - p)/pV^ , 

rotor thrust coefficient. 

2 

doublet strength, ft /sec. 

number of the collocation points. 

subscript indicating surface normal direction. 

2 

free stream pressure, Ib/ft . 

2 

static pressure at a point on upper surface, Ib/ft . 

o 

static pressure at a point at lower surface, Ib/ft . 
rotor radius, feet. 

radial distance from a vortex center, 
viscous core radius of a vortical element, 
lifting surface area. 

the axial, radial and the swirl velocity components of a 
vortex point. 

free stream velocity, ft/sec. 



LIST OF SYMBOLS (Cont'd) 


V 

n 

^imom 

w 

mn 


x,y,z 




a 


o 


"tpp 

r 

AD 

Ail 

AV 

At 

AC 

P 

y 


n 


o 


p 

a 

4 


normal component of the velocity. 

uniform induced velocity of rotor (momentum)^ ft /sec. 

normal velocity component at a collocation point due 
to all known effects^ ft/sec. 

components of the lifting surface coordinate system. 

the vertical distances of the two ends (A and B) of a 
vortex element from the surface, respectively, feet. 

the angle-of-attack of lifting surface with respect 
to free stream. 

rotor tip path plane angle-of-attack, positive aft, 
blade coning angle. 

turbulent kinematic viscosity of the vortex core. 

2 

circulation, ft /sec. 

o 

rotor blade bound circulation, ft /sec. 
incremental doublet strength, ft /sec. 
incremental lift airload, Ib/ft 

tangential velocity increment due to singularity. 

time increment for one computational step. 

2 

equals 2(p^. - P^)/pVq . 
advance ratio, V^/J2R. 

distance along the direction perpendicular to the vortex - 
axis direction. 

free stream density, Ib-sec^/ft^. 


rotor solidity. 


LIST OF SYMBOLS (Cont'd) 


disturbance potential distribution, 
potential at a point on the lower surface, 
potential at a point at the upper surface, 
blade azimuthal angle, 
rotor rotational speed, radians/sec. 



II. GENERAL FORMULATION 


The basic tools utilized to determine the influence of rotor airflow 
on empennage excitations are the UTRC Rotorcraft Wake Analysis (Ref. 7) 
and a nonlinear lifting surface analysis for unsteady three-dimensional 
flow. The wake analysis is used to determine the main rotor wake induced 
effects near the empennage surfaces and to provide the positions and 
strengths of the vortices located near the empennage surface. The detailed 
description of the technical approach, applications, and correlation re- 
sults for the Rotorcraft Wake Anlaysis are presented in Ref. 7. The paths 
and spacings of the vortices which would Interact with the empennage sur- 
faces, that is, those trailed off near the aft = 0®) portion of the 
rotor and near the forward (\p = 180°) portion of the rotor are shown sche- 
matically in Fig. 2. Once the position of the interacting wake is defined, 
the problem is reduced to the prediction of aerodynamic loads on lifting 
surfaces in the presence of concentrated vortices. The analytical formula- 
tion of this problem is described below. 


Statement of the Problem 

The objective of the analysis is the prediction of the unsteady airloads 
on low to moderate aspect ratio lifting surface generated by the close passage 
of the concentrated rotor vortices moving with the free stream. The lifting 
surface is assumed to be operating at low angles of attack, with no flow 
separation and the position and the strengths of the interacting vortices 
are known at all instants of time. The analysis considers the unsteady incom- 
pressible flow over a swept, three-dimensional wing surface of arbitrary shape 
at a given angle of attack. Lifting surface theory is used to predict, in 
time domain, the wing loading distributions due to the movement of the inter- 
acting vortices near the surface. The analysis includes the unsteady and 
suction effects of the interacting vortices, the unsteady effects of far main 
wake, and the unsteady effects of the shed vorticity behind the lifting sur- 
face. The interacting (viscous line) vortices are represented by finite 
segments of vortical elements during their interaction with the lifting sur- 
faces. 
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Basic Equations 


Consider the coordinate system of a wing as shown in Fig. 3. The x-y 
plane describes the zero mean angle-of-attack plane of the wing and the 
z-axis is directed away from the upper surface of the wing. In this coor- 
dinate system, the positions and the strengths of the interacting vortices 
are defined as a function of time. Since the wing loading varies with time 
due to passage of these vortices, any change in wing circulation results 
in shed vorticity. The boundary condition for the problem is that there 
is no flow normal to the boundary surface. 

The flow over the lifting surface is assumed to be potential; and can 
be described by an unknown disturbance velocity potential <f> (x,y,z) which 
varies with time and is zero at points far removed from the lifting surface. 
Define Vj^(x,y,z) as the instantaneous total normal velocity at any point on 
the lifting surface due to all the effects. The various components of v^ 
may include the induced effects of the disturbance potential at other points 
on the wing, the effects of shed and trailing vorticity behind the wing, 
the free velocity component due to angle of attack, and the velocity induced 
by any other known vorticity elements such as the interacting main rotor 
wake. The satisfaction of flow boundary condition results in 


= V , on surface S. (1) 

9n n 

Thus, the problem becomes the determination of a disturbance potential, <{>, 
as a function of time in a stepwise manner caused by an airfoil moving at 
a constant velocity through an unsteady fluid and shedding vorticity at the 
trailing edge. It can be shown that assuming inviscid, irrotational, and 
incompressible flow the disturbance potential <j) satisfies the Laplace 
equation (Ref. 8) 


V^(j) = 0 (2) 

at each instant of time. Thus the specific problem is the solution of the 
Laplace equation for the disturbance potential, <(), at each instant of time 
such that the time varying boundary condition described by Eq. (1) is sat- 
isfied . 
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Solution of the Laplace Equation 


It can be shown that the potential function induced by singularities 
such as a source, or doublet, will identically satisfy the Laplace equa- 
tion and will vanish at infinity. Therefore, the solution of the Laplace 
equation is one of finding a singularity distribution on the surface S 
that satisfies the normal boundary conditions in Eq. (1). Since a source 
distribution alone does not produce any resultant lift, a doublet distribu- 
tion is required. 

If D denotes a surface doublet distribution whose axis is everywhere 
along the outward normal to the local surface S, the potential induced at 
any point P is given by 



where r is the distance vector from doublet origin on S to point P, and n 
is the unit vector outward normal on the elemental local surface ds, and 
S is the total surface on which doublets are distributed. In general, the 
surface S consists of the upper and lower surfaces of the lifting surface 
and the wake surfaces. 

Substitution of Eq. (3) into Eq. (1) results in the following integral 
equation 


1_ 

4tt 


3n 




v 

n 


(4) 


The integration of the above equation determines the required doublet 
distribution D. This is carried out numerically in the analysis. 


The Boundary and the Trailing Edge Conditions 

The flow tangency boundary conditions (Eq. (1)) should be satisfied on 
the true wetted surface of the lifting surface as has been done in Ref. 3. 
This results in computations of the surface pressures on the upper and lower 
surfaces of the airfoil separately. But for this analysis, it is adequate 
to compute only the pressure differences between the upper and the lower sur- 
faces. This is sufficient since no flow-separation is expected and thickness 
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effects are expected to be small. As a result, the no-flow boundary conditions 
are satisfied only along the mean line of the airfoil, thus the thickness and 
the curvature effects for the computations of the airloads have been ignored. 
This reduces computational time required and thus makes the analysis more 
efficient. Therefore, in this analysis the doublet distribution on the wing 
represents the distribution of the potential difference between the upper 
and lower surfaces, or 


D(x,y,z) = (}i^(x,y,z) - <|>j^(x,y,z) 


(5) 


So far, in the development of the equations, wake trailing behind the 
wing has not been considered. It can be shown, by using the Helmholtz vor- 
ticity theorem along the trailing edge of the wing, that the doublet 
strengths of the wake trailing and shed elements can be expressed in terms 
of the doublet strengths imparted to those elements as they leave the wing. 
Thus, the elements on the wake do not introduce any new unknowns. In this 
analysis the wing wake geometry is assumed to be planar and consisting of 
the trailing elements extending to infinity and a finite number of shed 
elements, as shown in Fig. 4. 


Discussion of Numerical Method 

In order to integrate Eq. (4), the surface of a wing is assumed to be 
divided into a finite number of quadrilateral elements, as shown in Fig. 4. 
The value of the doublet strength D is assumed to be constant over each 
surface element and that the axis of the doublet is directed along the 
local surface normal at each control point which is at the centroid of 
each element. A numerical procedure is then utilized such that Eq. (4) 
is satisfied at a finite number of points corresponding to the centroids 
of the various elements representing the wing surface. 

It can be shown, by analogy with electromagnetic theory, that the flow 
induced by a doublet distribution of density D over a given area is the 
same as that due to a vortex of strength D around its boundary. Therefore, 
numerically the doublet distribution on the surface of the wing corresponds 
to a network of vortex elements as indicated in Fig. 4. The Biot-Savart 
law is utilized to obtain the velocity field induced by the network of vortex 
elements. Further discussion of this approach and specifically the discus- 
sion of the computer procedure by which this is accomplished is provided in 
a later section (see Section III) . 
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Determination of Pressure Distributions 


Once the doublet distribution has been obtained, Bernoulli's equation 
is utilized for the determination of wing pressure loading under the assump- 
tion of potential flow. The loading may be described by a pressure coefficient 
defined as follows 

ACp = (pj^ - p^)/ Y PVo^ ( 6 ) 


Considering a differential element of doublet, AD, which represents the 
elemental circulation, and the velocity, AV, just above and below the thin 
elemental sheet at the mean line (average velocity V) of length, dx, the 
circulation is 

AD = 2 AV dx, or 2AV = — K/) 


If pressure at a distance far from the airfoil is p^, Bernoulli's equation 
between the upper and lower surfaces (Ref. 8) becomes 


I p + Po = P -^ + I P (V + AV)^ + pu (8) 

and 

2 1 2 
1 V 1- p = p -7— + T P (V - AV) + p (9) 

— p o o ot / ^ 


Combining Eqs. (5), (6), (7), (8), and (9) results in 



( 10 ) 
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Modeling of Nonlinear Suction Loads 


Suction lift is the component of lift which results from the low 
pressure region within the vortices when they are in the close proximity 
to a lifting surface. The close interaction of a concentrated free vortex 
with a lifting surface can be separated into two mechanisms; one, the 
influence of the vortex induced velocity field, and second, the effect of 
the viscous core on the near field pressure distribution of the lifting 
surface. The first effect can be easily accounted for by including it 
during the determination of the potential flow field. For example, the 
doublet lattice method discussed in the earlier sections adequately in- 
cludes the influence of the vortex induced velocity field. The modeling 
of the second mechanism, the nonlinear suction lift, is of significance 
because available experimental measurements clearly indicate that large in- 
cremental suction peaks are generated during the close wing vortex inter- 
action and these suction peaks cannot be adequately handled by the conven- 
tional wing theories. The following section describes the analytical formu- 
lation of the suction load model incorporated in this analysis. 

In general, the relative velocity at any point on a real (viscous) 
vortex element moving with the free stream can be described by three com- 
ponents, these being the swirl component w, the radial component v, and the 
axial component u. For example, all three of these components exist in the 
core of a vortex during its formation near a blade tip (or near a leading 
edge of a highly swept wing) . During its travel from the blade tip to 
empennage surface, the vortex element goes through a roll up process where- 
by the radial and axial components of velocity become negligible and the 
vortex segment can then be represented by an element of a Rankine type vortex. 
For a Rankine type of vortex (with circulation F) the tangential velocity 
component is that of a potential line vortex outside of the viscous core 
region, 


w = r/2TTr, for r > r^. 


( 11 ) 


and the velocity component is that of a rotating rigid body inside the 
viscous core region 


w = Fr/2TTr for r < r^, (12) 

c 

so that all of the vorticity is confined to the viscous core. If we consider 
a Rankine vortex moving with a velocity, V^, in the vortex outer region, 
the radial pressure gradient of the vortex is solely balanced by the cen- 
trifugal force term. 
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(13) 


pw^/r = dp/dr 


The above equation is obtained from the Navier-Stokes equations described 
in cylindrical coordinates where the viscous and convective terms are 
neglected for r > r^. The integration of the above equation gives 


or 


p = p 


cx> 

- pr^/4ir^r^ dr 


(14) 


7 9 2 

p - p^'pr /Sir r 


(14) 


Defining the static pressure coefficient C as 

P 


2 

C = 2(p„ - p)/pV 
P o 

then for a vortex in a free stream 

C = (r/2irrV )^ 

P O 


(15) 


(16) 


The centerline of the vortex is assumed to be a streamline. Utilizing the 
method of images to satisfy the condition of no flow on the flat surface, 
the incremental ACp at any point (r) on the surface is given by 


or 


where 



9 2 

AC = A(r/V y/r 
P o 


A = 1/2tt 


(17) 


(18) 
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The above relationship corresponds to a simplified wing vortex inter- 
action model as shown in Fig. 5. For the general case where the vortex 
interaction with a wing is of an arbitrary nature as shown in Fig. 6, the 
incremental suction leading is obtained in an approximate way as discussed 
below. 


Defining AB, as a vortex element of projected length with a constant 
circulation T lying over a wing (see Fig. 6). Its distance from the wing 

(Z) varies linearly from point A to point B and can be expressed as 

The variation of suction loading in the n direction is such that the larger 

the distance Z, the more spread out the loading in the n**direction. In other 

words, the distance ti^ at which this loading goes to zero in the n^direction 
can be described by the approximate relationship 


n 

o 

The value of m varies between 1 and 2. 
along segment AB can be written as 


= mZ (20) 

The variation of suction load 


where 


AC = K/(n^ + z^) 
pv 


K = A(r/V^)^ 


( 21 ) 


( 22 ) 


The integrated load (Ail) on the wing due to segment AB is given by 


n 


dC 


(23) 


Evaluation of Eq. (23) results in 


-1 


A£ = K(2 tan m) £n ~ if \ ^ (24) 
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and 


A£ = K(2 tan /Z^ if Z_ = Z. 

O A D A 


( 25 ) 


I 

Similarly, the evaluation of the center of pressure of this load Ai. 
(denoted by point P in Fig. 6) is given by 


l/(£n|Zg/Z^|) - Z^/(Zg - Z^) if Zg Z^J , (26) 

where 

Equations (24) through (27) are used in this analysis to compute the suction 
loads. The value of constants A (Eq. (17)) and m (Eq. (20)) can be varied 
through input. Also, the values of or Z are always assumed to be greater 
than the viscous core radius of the vortex element AB. 
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III. DESCRIPTION OF COMPUTER PROGRAM 


The analysis described in Section II has been incorporated into the computer 
program RIEVA (Rotor Induced Empennage Vibration Airloads) . The flow chart of 
the key computational steps is shown in Fig. 7. The basic assumptions of the 
nonlinear lifting surface theory previously described are inherent in the 
program. 

(a) It is assumed that the strength and position of the interacting 
vortex elements are known at all instants. In general, the circulation and 
the coordinates of the interacting vortical elements are input to the program. 
The program does have a capability, however, to compute internally the 
coordinates of interacting wake elements corresponding to a classical undis- 
torted rotor wake. The equations for the description of a classical wake 

are given in APPENDIX A. 

(b) The distortions of the interacting vortex elements due to loading 
on empennage surfaces are not included. 

(c) The wake behind the stabilizer/fin trailing edges is assumed to 
consist of semi-infinite trailing elements and finite shed elements, all 
lying in the plane of the corresponding lifting surface (see Fig. 4). No 
roll up effects or distortions of this wake are included. 


Numerical Method of Solution 

The numerical procedure involved in applying the lifting surface method 
consists of dividing the surface into a number of appropriately shaped boxes. 
The number of these boxes is arbitrary and it is controlled through input. 

The magnitude of the doublet strength D over each box is assumed to be uni- 
form. The total velocity induced perpendicular to the surface at a colloca- 
tion point (centroid of the box) consists of that due to the vorticity of 
all other boxes on the surface, the effects of the concentrated interacting 
vortices, the velocity induced by the far rotor wake and that due to all the 
trailing and shed vorticity elements starting at the trailing edge of the 
lifting surface. When the flow tangency requirements on the surface are 
satisfied at each time step, the problem of calculating the doublet strengths 
is reduced to one of solving N linear algebraic equations, where N is the 
number of boxes on the lifting surface. Specifically, if aero- 

dynamic influence coefficient at the centroid of the box mn du4 to the effect 
of box ij and w is the normal component of the velocity at collocation 
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point mn due to all known vorticity elements and due to angle of attack, the 
satisfaction of the boundary conditions results in the following relationship 


Z Z 

A . .D. . = w 
^ j mnij ij mn 


( 28 ) 


Since the coordinate system chosen in this analysis is fixed to the wing 
at all times, the influence coefficients matrix is computed only once 

and stored for use at subsequent steps. 


Discussion of Options 

One of the essential features of the RIEVA analysis is that the total 
main rotor wake has been split into two parts for the computations of in- 
duced velocities at the collocation points. One is the interacting elements 
represented by a small number of concentrated vortices which come very close 
to empennage surfaces causing high harmonic induced effects. The other is 
the induced velocities due to the rest of the wake (far wake) varying at a 
relatively low magnitude. As a result, far wake induced velocities are com- 
puted at relatively large intervals of time and their magnitude at intermediate 
intervals is calculated by interpolation. As mentioned earlier, this pro- 
cedure results in savings of computational time without the elimination of 
high frequency effects. Some additional advantages are discussed below. 

As indicated by the Biot-Savart law, the magnitudes of the high frequency 
airloads induced on the empennage surface depend upon the relative distance 
between the surface and the interacting vortex elements (besxdes the circula 
tions of the vortex elements). For a helicopter rotor in forward flight, 
only the undistorted wake can be described by analytical means. It is very 
difficult to compute accurately the distortions of the wake elements while 
they move from rotor blade to the empennage surface. In fact, some of these 
interacting elements go through distortions of very large magnitude due to 
encounters with rotor hub, fuselage, etc. Therefore, the positions of the 
interacting vortices with respect to the empennage surfaces can only be 
approximated. In the RIEVA analysis, since the interacting rotor wake is 
represented by a limited number of vortices (normally blade tip vortices), 
it becomes relatively easy to modify their geometry and thus compute the air- 
loads corresponding to the modified wake geometry. As a result, once the 
approximate position of the wake is defined, the interactions could be real- 
istically assumed to occur anywhere within the envelope described by varying 
the distance between the empennage surfaces and the vortices. The analysis 
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is flexible enough to consider many types of interactions. Specifically, the 
interacting vortices can be arbitrarily moved closer to the empennage surface 
to determine the magnitude of the critical airloads corresponding to the 
extreme cases of wake distortions. Thus, wake distortion effects can be 
accounted for indirectly and thereby the analytical computations of the 
distorted rotor wake geometry are avoided. 

The various options available in the RIEVA analysis for describing the 
interacting wake coordinates and the rotor induced velocities are described 
below. 

For interacting wake geometry description, the user has the option to 
input the coordinates of interacting wake elements at each instant of time. 
This is particularly useful for those cases where the distorted wake geo- 
metry is known. Alternately, the user can utilize the wake coordinates 
computed internally by the RIEVA analysis. For a given flight condition 
and rotor parameters, the program computes the wake coordinates corresponding 
to the classical undistorted skewed helical wake. 

Also, for the inclusion of far wake effects, there are two options 
available. For one option the user can compute the induced velocities at 
empennage points due to all or part of the rotor wake by utilizing a 
Rotorcraft Wake Analysis (Ref. 7) and input these velocities to the RIEVA 
analysis. Alternatively it may be assumed that the far wake induces a 
uniform velocity at the empennage points and its magnitude is described 
through input. This approximation may be valid for some high speed flight 
conditions. The coupling of program RIEVA with the Rotorcraft Wake Analysis 
is further discussed in the following section. 


Coupling With Program F389 

As indicated in the flow chart (Fig. 7), the Rotorcraft Wake Analysis 
(Program F389, Ref. 7) may be utilized to obtain some of the inputs required 
for the execution of program RIEVA. These inputs may be divided into three 
convenient groups as described below. 

1. The circulations of the interacting vortex elements can be deter- 
mined by utilizing the Program F389. For a given rotor configuration, with 
known control inputs. Program F389 computes the radial and azimuthal distribu- 
tion of the blade circulation, r^(r,i|^). Once the distribution, F^(r,i^), is 
prescribed, the circulations of the interacting vortices can be very easily 
estimated by using a simple roll up theory. Normally these interacting vor- 
tices will consist of only the blade tip vortices that are trailed off the 
fore and the aft portions of the rotor disk. 
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2. The induced velocities at the empennage collocation points due to 
main rotor wake may be computed by utilizing the Program F389. Normally, 
these induced velocities are computed at coarse intervals of time and the 
effects of the interacting elements (elements in group 1. above) are ex- 
cluded from these computations (this is done to avoid redundancy) . 

3* If desired, the coordinates of the interacting elements can also 
be obtained from the output of the computer program F389. 


Interacting Vortex Geometry Modifications 

This section describes the procedures utilized by the RIEVA program to 
automatically modify the geometry of the interacting wake elements in case 
these elements intersect the empennage surfaces • This is carried out in a 
very approximate manner. The use of this procedure is optional and it is 
specified through input to the program. 

Figure 8 illustrates the modifications involved. The modifications are 
confined to junction points (where the various elements of vortex meet) . The 
junction point nearest to the wing surface (only if the vortex element is cut) 
is modified, as shown in Fig. 8. The junction point is always kept at least 
one core radius (viscous core radius r^) away from the surface. Also the 
modified junction points at locations above the wing are not allowed to 
penetrate the surface at subsequent time steps, but instead they are allowed 
to slide over the wing surface (or even some distance beyond the trailing 
edge, this distance being specified through input). 

The modifications described above involve the changing of only the 
Z-coordinate of the affected vortex elements. No attempt is made to deter- 
mine the vortex displacements from their original position due to the wing 
vorticity field. 


Future Capabilities 

Some of the improvements that may be carried out in the future to enhance 
the capabilities of the RIEVA program are discussed in the following section. 

1. At present a considerable amount of user judgement is required to 
divide the main rotor wake into two parts; one, the interacting elements and 
two, the far wake elements. A procedure should be developed and then in- 
corporated in the RIEVA analysis whereby this splitting of wake is carried 
out automatically within the program. 
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2. The nature of the interaction between the main rotor vortices and 
the vertical tail is such that significantly large streamline velocities 

may be induced on the tail surface by these interacting vortices. This span- 
wise variation of the strearawise velocity results in a significant amount 
of vibratory airloads being generated by the vertical tail. At present, 

RIEVA analysis does not account for these effects. 

3. The coding and the dimensions of the RIEVA program should be modified 
so that the airloads acting on the vertical fin and the horizontal stabilizer 
can be computed simultaneously. 

4. Flow separation effects should be included in the RIEVA analysis. 
Normally, when a vortex is close to the surface, flow separation occurs on 
the upwash side of the vortex. 
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IV. PROGRAM RIEVA DOCUMENTATION 


The computer program for the determination of rotor induced empennage 
vibration airloads (RIEVA) incorporates a number of key options. The 
inclusion of these options results in computer efficiency and provides 
flexibility on the type and source of the interacting wake elements. 

The program documentation includes 

(i) a brief description of program subroutines, 

(ii) a detailed description of input data, 

(iii) a brief description of program output. 


Subroutine Descriptions 

A brief description of each of the subroutines is given below in 
alphabetical order. 

Subroutine COLLQC 


This subroutine computes the aerodynamic influence coefficient matrix 
A at the first time step and stores it on unit 10 for use at subsequent 
tTSe*^ steps. It also computes the coordinates of the collocation points and 
prints out the coordinates, if requested. 

Subroutine DMAT 

This subroutine computes the solution of a set of simultaneous algebraic 
equations of order N. 

Subroutine GETVIS 


This subroutine interpolates input induced velocity components corre- 
sponding to the far rotor wake at the intermediate time steps. 

Subroutine HARMPM 


This subroutine computes the harmonic coefficients of the wing 
aerodynamics loads. 
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Sub rout ine LOADER 


This subroutine reads in the loader input data. 

Subroutine LOADS 

This subroutine computes and then prints out the wing aerodynamic loading 
distribution AC at each time step. Also, the total wing aerodynamic loads 
are computed an8 printed out. 

Subroutine MAIN 

This is the main routine. It calls subroutines LOADER, PRNT, WLINPT, 
COLLOC, VRTXG, SLOAD, GETVIS, VELKNW, SOLVEG, LOADS and HARMPM and prints 
out the harmonics of the wing airloads. It also prints out the instantaneous 
geometry of the interacting vortices. 

Subroutine OUTPUT 


This is a general subroutine which prints out the chordwise and spanwise 

distribution of the various quantities (e.g., doublet distribution D, pressure 

distribution Ac' , etc.) 

P 

Subroutine PRNT 


This subroutine prints out all the input loader data. 
Subroutine SLOAD 


This subroutine computes the nonlinear suction airloads, when requested. 
Subroutine SOLVEG 

This subroutine sets up the simultaneous set of algebraic equations to 
be solved at each time step, and also, outputs the computed doublet distribu- 
tion. 

Subroutine UVW 


This subroutine computes the induced velocities due to a vortical element. 
Subroutine VELKNW 


This subroutine computes the induced velocities at the collocation points 
due to the interacting vortices and also due to the wing wake trailing and 
shed elements. 
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Subroutine VRTXG 


This subroutine computes the coordinates of the interacting main rotor 
vortices corresponding to a classical undistorted wake model. 

Subroutine WLINPT 


This subroutine reads in the main rotor far wake induced velocities and 
converts them to the wing coordinate system. It also prints out the 
distribution of the normal component of these far wake induced velocities 
corresponding to each time step. 

Input Description 

The required input to the program consists of the following major punched 
card data blocks, in order of loading; 

I. Wing (Empennage), Vortex Geometry and other Solution Control 
Data (in Loader Format) 

II. Induced Velocity Data due to Far Main Rotor Wake (Optional) 

III. Interacting Vortex Geometry Data (Optional) 

Details for preparing these blocks of data are given in the following 
sections. 

I. Empennage, Vortex Geometry and Other Solution Control Data 

This block of data is split into the following three groups. 

(a) Solution Control Parameters 

(b) Empennage 

(c) Vortex Geometry Parameters (Optional) 

The group (a) contains the data that determines which of the various 
options are to be exercised. The data in this group are loaded in an 
array (SCP). The group (b) contains data that describes the empennage sur- 
face geometry and the free stream parameters. This group of data is 
stored in WCP array. The group (c) contains the rotor data needed to des- 
cribe the interacting main rotor wake geometry based on a classical undistorted 
model of the rotor wake. This group of data is stored in array VCP. If 
internal computation of the interacting wake geometry is not required, this 
group of data may be omitted. 
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The order in which the groups of data are loaded is arbitrary. The 
format for these data is as follows; 

A N L DATA(L) DATA(L+1) .... DATA(L+4) (Al, II, 14, 5F12.0) 

where A, in Column 1, represents the letter S, W or V corresponding to arrays 
SCP, WCP or VCP. N is the number of data items to be input on the card 
Column 2; N must not exceed 5. L is the location or identifying number of 
the first data item on the card Columns 3-6, right adjusted. DATA(LH-j) 
represents the various data items on the card, Columns 7—18, 19—30, 31—42, 
43-54, and 55-66, in floating point format. The locations for the various 
data are listed below along with definitions and other pertinent comments. 
Note that some locations are intentionally left blank. 


(a) ^CP Ar^a^ 
Location Item 


1 


N, 


W 


2 MYES 

3 NSHED 

4 MSYM 

5 NVORT 

6 IVGEOM 


Description 

The number of lifting surfaces. For 
horizontal stabilizer N^=2., one for 
each semi-span. For vertical tail, 

Ny=l. At present both stabilizer and 
fin cannot be modeled simultaneously. 

MYES-1., doublet lattice method is 
applied. MYES=0-, no doublet solution 
required; this option is used if only 
nonlinear suction loads are desired. 

Specifies the number of shed elements 
to be included in the wing wake. 

Range is 0. to 5. (see Fig. 4). 

Normally MSYM=0. When MSYM=1.> the pro- 
gram assumes symmetrical loading on both 
the right and left hand sides of the 
stabilizer. 

Specifies the number of interacting 
vortices. Range is 0. to 8. 

Option for the description of the inter- 
acting wake elements. IVGEOM=0, rotor 
wake data specified via loader inputs 
in VCP array (group c) . IVEG0M=1, 
interacting wake data input at each time 
step. This is described separately (data 
block II, Page 29). 
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Location 


Item 


Description 


7 NFARF 


8 MODVG 


9 EXT 


10 INSUC 


11 NT 


12 

AT 2 

13 

MHDO 


Option for including the far wake 
induced velocity effects. If NFARF=1, 
this component of the induced velocities 
is to be input as described separately 
(data block II). If NFARF=0, no data 
are read in. If desired, under this 
option (NRARF=0) , these effects may 
be accounted via input in location 19 
of SCP array. 

Option to modify the interacting vortex 
geometry if the vortex gets cut by the 
empennage surface (see Section II, 

Page 13), For MODVG=0, no modifications 
are carried out. 

If M0DVG=1 (previous location) , the 
area of the wing surface over which the 
modifications are to be carried out is 
extended via input in this location, i.e. , 
the modifications of the vortex are carried 
out even beyond the wing surface. EXT 
is nondimensionalized with respect to the 
mean chord length of the wing. 

Option to include the computations of 
nonlinear suction loads. INSUC=1, loads 
are computed. INSUC=0, no loads are 
computed . 

Specifies the total number of time steps 
for which computations are to be carried 
out • 

Specifies time increment per step, secs. 

Option to perform an harmonic analysis 
of the computed air loads over one period - 
For MHD0=1, Yes and for MHDO=0, No. The 
number of points per period is specified 
in location 15. 


24 



Location 


Item 


Description 


14 NH 


If MHD0=1 (location 13), the input spec- 
ifies the number of harmonics desired. 


15 

16 


17 


18 


19 


20 

21 


22 


23 


NT2 

A 


T 

n 


e 


fm 


F 

V 


LPCOLL 


LPVRTX 


LPINDV 


Specifies the number of steps per period. 

The constant A in the suction airload 
model (Equation 18) if different from 
1/27 t^ (default value) . 

The parameter (2 tan ^ m) in the suction 
airload model (Equation 25) if different 
from 2.21429 (default value). 

Factor to modify the first shed element 

from the TE. Included for possible future 

use (Reference 10). A nonzero positive 

value of e^j^ (0. to 1.) models the first 

shed element closer to the TE by the 

factor e^ . 

fm 

If the effects of the far rotor wake 
are represented by a constant uniform 
induced velocity over the wing, its 
magnitude is input here. It is non- 
dimensional with respect to 
(positive downwash) which should be 
input in location 3 of VCP array. 

Not used. 

Print option for coordinates of the 
collocation points over wing. LPC0LL=1, 
Yes. LPCOLL=0, No. 

Option to print vortex geometry 
coordinates at each time step. LPVRTX=1, 
Yes. LPVRTX=0, No. 

Option to print the total (near plus 
far wake) induced velocity distribution 
at each time step. LPINDV=1, Yes. 
LPINDV=0, No. 


25 



Location 


Item 


24 LPDBLT 


25 LPPRCP 


26 LPFARW 


27 LPSUCD 


28 IDEBUG 


Description 


Option to print the doublet distribution 
at each time step. LPDBLT=1, Yes. 
LPDBLT=0, No. 

Option to print the airload distribution 
AC at each time step. LPPRCP=1, Yes. 
LpPrCP=0, No. 

Option to print the inputted (data block 
II) far-wake distribution. LPFARW“1, 
Yes. LPFARW=0, No. 

Option to print the computed nonlinear 
suction airloads. LPSUCD=1, Yes. 
LPSUCD=0, No. 

Option to print debug parameters. 

Always input IDEBUG=0. 


(b) WC^ 

1-2 Xj^,Y^ 


3-4 

5-6 

7-8 

9 


X2.Y2 


X3»Y3 


^ 4 >\ 


N 


10 


N 


11-20 


(") 


X and Y coordinates of corner 1 of the 
first surface in the wing coordinate 
system illustrated by Figure 4 (feet) . 

X and Y coordinates of corner 2 (feet), 

X and Y coordinates of corner 3 (feet). 

X and Y coordinates of corner 4 (feet). 

Number of chordwise stations for the 

lattice of the first wing. 

Number of spanwise stations for the 
lattice of the first wing. 

Specifies the parameters for the 
second wing in the same order as for 
the first wing (in locations 1 through 
10 ). 
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Location 


Description 


Item 

21 V 

o 

22 

23 n 

24 P 
(c) VCjP 


Magnitude of free stream velocity, ft/sec. 
The angle of attack of wing, deg. 

Mean reference chord, ft. 

2 A 

Density of free stream medium, lb-sec /ft^. 


Input of this data is optional. Input only if the value in location 6 
of SCP array is 0, and the value in location 5 of SCP array is greater than 

0 . 


Location Item 


Description 


1 


Blade tip speed, ft/sec. 


2 R 


Blade radius, ft. 


3 


4 


5 


6-8 


9 


10 


V 

imom 


“tpp 


8 

o 




IFB 


Uniform momentum induced velocity, 
ft/sec. (positive for inflow or 
positive thrust) 

Tip path plane angle of attack, deg. 
(positive aft) . 

Rotor cone angle, deg. 

X, Y and Z coordinates of rotor hub, 
specified in the wing coordinate 
system, ft. (see Figure 3). 

Eddy viscosity of the vortex core, 
ft^/sec. Required to compute the 
increase in core radius with time 
(see Equation B-1 , in Appendix B) . 

Option for the inclusion of blade 
bound circulation effects by modelling 
it as one of the elements interacting 
with empennage. If the blade passes 
very close to the empennage, IFB should 
be input as 1. 
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Location 


Item 




12 


N. 


13 


14-20 

21 


22 

a 

23 AiJ; 

24 Ai|^ 

ru 


25 r 

CO 


26 


27 A^2 


28-29 


Description 

Nondimensional factor (between 0, to 
1.) to multiply the maximum blade cir- 
culation in order to get effective circula- 
tion. 

Number of rotor blades. 

This is the input required to make the 
time reference of data block II compatible 
with the time reference of the interacting 
wake geometry. is the reference step 

number in data block II that corresponds 
to starting time (t=0) of wake geometry 
computations . 

Not used. 

Azimuthal location at which the first 
point of the vortex is trailed, deg. 

(see Equation A-4) . 

Age of the first point of the vortex, deg. 
(see Equation A-5) . 

Azimuthal increment, deg, (see Equation 
A-4) . 


Time for the vortex to roll up to the 
initial core radius r^^ specified in 
location 25, deg. 

Initial core radius of the vortex, ft. 

Radial location on blade from which 
the vortex is trailed (nondimensional) . 
r^=l corresponds to tip vortex (see 
Equations A-1, A-2) , 

Azimuthal increment used in the compu- 
tation of the far wake effects, deg. 
Required for the aft forming vortex only. 

Not used. 
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II, 


Location 

Item 

Description 

30 

N 

V 

Number of segments representing the 
vortex. Maximum of 9. 

31-39 

r 

V 

Circulations of segments, ft^/sec. 

40 

- 

Not used. 

41-59 

(") 

Corresponding data for the second 
vortex similar to those for the 
first vortex in locations 21-39. 

60-179 

(") 

Similar data for the remaining six 
vortices . 

Induced Velocity Data (Optional) 


This block of data is required only if the value of location 7 of the 
SCP array (see data group a) is 1,0. F389 (Ref, 7) (see Section III, Page 

17) generates this data in a format which is compatible with the input 
requirements of RIEVA program. This block of data consists of three com- 
ponents of the instantaneous induced velocity at each of the collocation 
points. These induced velocities are described in the rotor tip path plane 
coordinate system (see Appendix A), Also, these induced velocities are 
described for various time steps; the total time covered being exactly 
the time for one blade passage. The first card of this data block contains 
three quantities; LS, AT, and NTOT is the format (12, F20.0, 15). LS is 
the total number of time steps, AT, is the time in seconds corresponding 
to each step and NTOT is the total number of collocation points. The 
total number of quantities (N^) contained in the data (following the first 
card) is 

N , = 3 X LS X NTOT 
a 

where 3 represents the three components of induced velocity. 

These N^j numbers are input in order with six quantities per data card, 
in the format (6F12.0). 

Ill • Interacting Vortex Geometry Data (Optional) 

This block of data is required only if the value of location 6 of 
SCP array (see data group a) is 1.0. This input data block describes 
instantaneous locations and strengths of all the interacting vortex 
elements. This option may be utilized if the Interaction between the 
empennage and the vortices is of non-periodic (transient) nature. Also, 
if the main rotor wake distortions are known, the wake data for these 
distorted elements may be input via this block. 


For each of the vortices, the following quantities are input. 
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Format 


Card // 


1 

2 


3 


Quantity or Quantities 

(Number of vortex points) (15) 

x,y,z,R,Gv (for first point) (5F10.6) 

where: x,y,z are the coordinates 

in the wing coordinate system, 

R,Gv are the viscous core radius 
and the circulation respectively, 

x,y,z,R,Gv (for second point) (5F10,6) 


N^+1 x,y,z,R,Gv (for N-^th point) (5F10.6) 


The above information is input for all the NVORT (see location 5 
of SCP array) vortices. 

The above information is repeated for each time step until data for 
all time steps has been input. 


Output Description 

A brief description of the printed output generated by the RIEVA 
program is given here. The amount and the nature of the output for a 
particular run is controlled through the exercise of various print 
options specified via input in locations 21-27 of the SCP array (see 
Section IV, Page 25) . 

The RIEVA program output can be classified into the following ten 
major categories (in order of output). 

1. Listing of Input Loader Data 

2. Listing of Input Far-Wake Induced Velocities 

3. Coordinates of Collocation Points 

4. Instantaneous Description of Interacting Vortices 

5. Instantaneous Induced Velocity Distribution 

6. Instantaneous Doublet Distribution 

7. Instantaneous Airload Distribution 

8. Instantaneous Integrated Airloads 

9. Nonlinear Suction Loads (if present) 

10. Harmonics of Airloads 
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The categories 2,, 3., 5., 6., and 7* listed above involve the 
chordwise and spanwise variation of the various quantities and are 
printed out in a similar format by the subroutine OUTPUT. Category 3. 
provides the coordinates of the various collocation points over all 
the surfaces, and the numbers in the output categories 2., 5., 6., 
and 7. represent the values of the various quantities at the corres- 
ponding points (listed in category 3.). 

Most of the output of RIEVA analysis is self-explanatory; only a 
brief description of each category is provided. No samples of the 
output are included here. Also it should be noted that categories 4. 
through 8. are output at each time step. 

1. Listing of Input Loader Data 

This is output by Subroutine PRNT, which lists the values of all 
the input loader data (discussed in Section IV. under Input Descrip- 
tion: Data Block I.). 

2. Listing of Input Far-Wake Induced Velocities 

This is output by Subroutine WLINPT, which lists the input induced 
velocities at the empennage points. First the total number of time 
intervals and the time (in seconds) corresponding to each interval are 
printed out. 

3. Coordinates of Collocation Points 


This output by Subroutine COLLOC, which lists first the x-coordinate 
of the control points and then the y-coordinate of the control points in 
the wing coordinate system. Before these collocation points are listed, 
the coordinates of the corners of the various control surfaces modeled 
(along with some flight parameters) are given first. 

4. Instantaneous Vortex Geometry 

This is output from program MAIN and consists of the description of 
position and strength of the various elements of the interacting vortices. 
First the step number and the time (in seconds) corresponding to that 
step are printed out. Next, for each junction point of the interacting 
vortices, x, y and z coordinates (in wing coordinate system) are listed. 
Also the circulation (ft**2/sec.) and the viscous core radius (feet) of 
each of the elements are listed. Finally a quantity ZMOD is listed. If 



the value of ZMOD is 1,0, it indicates that the z-coordinate of the junction 
corresponding to that point has been altered according to the procedure 
given in Section III, Page 18). 

5 . Induced Velocity Distribution 

This is output by Subroutine SOLVEG, which describes the instantaneous 
total induced velocity (due to far wake plus interacting wake) distribution 
(positive downwash) . 

6 . Doublet Distribution (listed as 'SOLUTION GAMMAS^) 

This is output by Subroutine SOLVEG, which lists the computed doublet 
distribution over the control surfaces, A positive value corresponds to 
a clockwise direction of the vortical elements around the box (see Fig. 4.). 

7. Attached Flow Airload Distribution 


This is output by Subroutine LOADS, which describes the unsteady AC^ 

(ACt^=C ' - C ) or airload distribution over the surface. 

F pj^ pu' 

8. Integrated Airloads 


This is output by Subroutine LOADS. At the end of the output of each 
step, the integrated loads for the total surface are listed. Both, the 
nondimensional coefficients and C , and the dimensional loads L 

(lbs), PM (lb. -ft.) and RM (lb. -ft.) are printed out. 


9. Nonlinear Suction Loads (if present) 

This is output from program MAIN and consists of the integrated 
suction loads (AC^y, AC^^^ and AC^^) for all the time steps over one 
blade passage. These are listed under the heading "Suction Loads 
Variation". 


10 . Harmonics of Airloads 


This is printed out from program MAIN. The airload (C^> or Cj^) 
variation over one blade passage is converted into mN^ per rev (N^ ~ 
number of blades, m=l,2,3,etc) harmonic airloads. For example total 
is expressed as 
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M 

*= a + Z a cos (mN, + b sin (mN i(^) 
L o m b b 

m=l 


or 

M 

C « a + Z r sin (mN, + (() ) 
L o tn ^ b^ 

in=l 


For each load the quantities a^, [a^, ^ ^ m=l, . . . ,M) ] are 

printed out. Maximum value of M is 4. The^various loads that are harmon- 
ically analyzed are listed below in the order of their print out: 


!• Ct , and Cr. for control surface 1. 

L M ^ 

2. C * C and for control surface 2. 

L M R 

3. Cj , C and C for total (1+2) surface 

^ M R 

4. AC^, and AC representing the incremental airloads due 

tv MV RV 

to nonlinear suction. 
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V. RESULTS AND DISCUSSION 


This section presents and analyzes the results obtained from the 
RIEVA (Rotor Induced Empennage Vibration Airloads) program for some 
realistic check cases, 

The basic check case chosen involved the computation of the unsteady 
airloads acting on the Black Hawk horizontal stabilizer surface for a high 
speed forward flight condition of 172 knots. Calculation of the strength 
of the interacting vortex elements and the induced velocities at the 
stabilizer points due to the rotor far wake were necessary before RIEVA 
could be executed. An undistorted classical wake model and trimmed for- 
ward flight were assumed for the calculations. The various flight and 
control parameters for this case are listed in Table I. The variation 
of the maximum blade circulation with blade azimuth position W 

is shown in Fig. 9. The total rotor wake has been split into two parts, 
one, the interacting wake elements and two, the far wake as shown in 
Fig. 10. The interacting wake elements comprised the segments of the blade 
tip vortices trailed off the forward and aft portions of the rotor (30 deg. 
on each side of the centerline or the x-axis) . The circulations of these 
tip vortex elements were the values of in Fig. 9. 

The far wake induced velocities at the stabilizer points were computed 
(using Rotorcraft Wake Analysis F389) at time intervals corresponding to Aip 
(^ * At) of 15 degrees. For one representative point on the stabilizer 
surface (at quarter chord point on the centerline) , the variation of the 
far wake induced velocity with time is shown in Fig. 11. Also shown in 
Fig. 11 is the variation of the total induced velocity (due to interacting 
wake plus far wake) with time. Examination of Fig. 11 reveals that most 
of the high frequency variation of the total induced velocity is caused by 
the interacting wake elements. This result substantiates the assumption 
that the computations of the far-wake induced velocities are only necessary 
at coarse intervals of time. 

RIEVA was used to calculate the instantaneous positions of the inter- 
acting vortex elements at fine intervals of time. For the basic case the 
total time corresponding to one blade passage (T = 0.05814 seconds) was 
divided into 12 equal steps and the computations were carried out every 
.004845 seconds. At three such instants (steps 1, 7, and 12) the paths 
and spacings of the interacting rotor blade tip vortices with respect to the 
horizontal stabilizer are described in Figs. 12, 13, and 14, respectively. 

In these figures, the vortices numbered 1, 2, and 3 are those trailed off 
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the forward portion of the rotor (i|j = 210 to 150 degrees) while the vortices 
numbered 4 and 5 are those trailed off the aft (\(; = 30 degrees to -30 degrees) 
portion of the rotor. The integrated unsteady airloads of the stabilizer 
(at zero angle of attack) are shown in Fig. 15 as a function of time (i(^=f2t) . 

In this figure C ^ is the lift coefficient (normalized to x wing 

area). is the pitching moment coefficient about the quarter chord 

(positive up), and C^rj. is the rolling moment coefficient about the stabilizer 
centerline (positive for right side going down) ; both of these moment co- 
efficients are normalized (.5pV^2 x wing area x mean chord). As expected, 
the airloads C^.^, and C ^ are periodic with each blade passage (or 
every 90 degrees for four-bladed rotor) . The computations have to be carried 
out for more than one blade passage (as shown in Fig. 15), since the 
strengths of the shed elements were not known at the start of the computa- 
tions. Once the shed vorticity develops (after 5 to 6 steps), the results 
indicate excellent periodicity with each blade passage (as indicated by 
the results shown in Fig. 15). The airloads’ variation over one period was 
converted into the 4m per rev (m = 1, 2, 3) harmonic airloads and result 
in the following vibratory airloads on the stabilizer: 


Clt = - 0.184 + 0.127 sin (4i(; - 90*^) 
+ 0.068 sin (8i|; + 168^) 
+ 0.032 sin (12ij> - 36®) 


= - 0.0032 + 0.0187 sin 
MX 

+ 0.0161 sin 
+ 0.0112 sin 


(4i^ + 129®) 

- 8.3®) 
(12i{j + 144®) 


C 


RT 


0.0058 + 0.0773 sin (4ip + 195®) 
+ 0.0750 sin (8iJ; + 98®) 
+ 0.0400 sin (12i|; - 31®) 


The airloads described above are the total loads being applied to the 
stabilizer of the helicopter and these can be used as input to the coupled 
rotor/airframe (Ref. 6) vibration analysis. The computer program RIEVA 
also provides a detailed description of these airloads (chordwise and span- 
wise variation of AC ) . These airloads may be utilized to study the aero- 
elastic vibrations o¥ the empennage itself. As an example, for this check 
case some calculated airload distributions are shown in Figs. 16 and 17 (for 
three spanwise stations; at 6% and 56% stations on the right stabilizer and 
at 56% station on the left stabilizer). Figure 16 describes these airload 
distributions at the time step No. 1, the instantaneous vortex geometry for 
which is illustrated in Fig. 12. The chordwise variation of the airload 
at this instant resembles the conventional airload distribution observed on 
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an airfoil at a negative angle of attack. But when a concentrated vortex 
is directly above the stabilizer (such as vortex number 2 in Fig. 13), 
these chordwise airload distributions look totally different. This is 
clearly indicated by the distribution of the airloads described in Fig. 17, 
which correspond to the time for which the instantaneous vortex geometry is 
described in Fig. 13. These deviations are obviously due to the unsteady 
effects caused by the presence of the moving vortex. 

As indicated by Figs. 12, 13 and 14, the closest distance between 
the stabilizer surface and the interacting vortices is approximately 1.5 
feet. If this distance were smaller (say less than three times the 
viscous core radius of the vortex) , either due to wake distortion or due 
to change in tip path plane angle of attack, an additional component of 
the airloads results. As discussed before, this component of the air- 
loads is due to the nonlinear suction effects that occur when the vortex 
passes very close to a surface. 

An additional check case was run to verify the modeling of the non- 
linear suction airloads in the RIEVA analysis. This case corresponds to 
the same flight conditions as described earlier (see Table I) except that 
the tip path plane angle was arbitrarily increased to -14.42 degrees. 

This was done to let the vortices off the front portion of the rotor pass 
very close (within 2 inches) to the horizontal stabilizer surface (see 
vortex #2 in Fig. 18). As a result of close stabilizer— vortex interaction, 
^^S^ificant nonlinear suction airloads result and as shown in Fig. 19 for 
one blade passage. These incremental suction airloads (AC , AC-^, AC ) 
have to be added appropriately to the potential flow loadable , Cj^’ an^Cj, 
shown in Fig. 20) in order to determine the total loads C , and C 
These total loads are shown in Fig. 21 for one blade passage. It should 
be remembered that the loads described in Fig. 21 have been computed under 
the assumption that no flow separation of any kind occurs. Normally flow 
separation occurs on the upwash side of the vortex during the close airfoil- 
vortex interaction. At present the RIEVA analysis does not account for 
these flow separation effects. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


A computer program (RIEVA) has been developed that predicts the 
unsteady aerodynamic forces that are imposed on the empennage surface 
due to its interactions with the main rotor wake. The program was 
exercised to determine the vibration airloads acting on the Black 
Hawk horizontal stabilizer for a high speed (y = 0.4) forward flight 
condition. The following are specific conclusions obtained from 
analysis of the Black Hawk simulation. 

1. The results demonstrate that it is possible to compute the 
high frequency empennage vibration airloads (up to 16 revolu- 
tion) efficiently (reasonable amount of computer time) by 
utilizing this analysis. This check case required only 4 to 
5 minutes computer time on an IBM 370 system. 

2. No numerical problems were encountered, and the airload 
harmonics converged quickly. 

The following future activities are recommended: 

* 

1. A correlation study should be carried out to compare the RIEVA 
predicted empennage airloads with test data. 

2, The RIEVA analysis should be extended to account for the span- 
wise variation in streamwise velocity (nonuniform free stream 
velocity) . The inclusion of this feature will result in the 
analytically correct computation of the vertical tail airloads 
and provide the capability of studying the blade finite tip- 
vortex interaction problem. 

* 

The draft copy of the present report was prepared in 1981. Since then a 
limited correlation study has been carried out wherein the RIEVA predicted 
stabilizer airloads have been compared with the CH-53A flight test data. 
The results indicate a reasonably good correlation between analytical 
predictions and the flight test data. These results are presented in the 
following Reference; 

Gangwani S. T. : Determination of Rotor Wake Induced Empennage Airloads. 

Presented at the American Helicopter Society National Specialists' Meeting 
on Helicopter Vibration, Hartford, CT., November 2-4, 1981. 
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Table I 


BLACK HAWK CHECK CASE PARAMETERS 


Free Stream Velocity, Knots 172.0 
Advance Ratio, y 0.39 
Blade Radius, Feet 26.83 
Thrust, Pounds 17,203 
C^/a .0902 
Tip Path Plane Angle, Degree -10.8 
Cone Angle, Degree 3.66 
Uniform Momentum Induced Velocity, FPS 6.67 
Rotor Solidity, a 0.0821 
Number of Blades 4.0 
Hub X-coordinate , Feet 28.25 
Hub Y-coordinate , Feet 0.0 
Hub Z-coordinate , Feet 6.25 
Area of Stabilizer, Square Feet 44.45 
Mean Chord of Stabilizer, Feet 3.083 
Stabilizer Angle of Attack, Degree 0.000 



APPENDIX A 


The wake equations for a rotor in forward flight, based on the 
classical undistorted helical wake model, are presented below. These 
equations are utilized in the analysis (RIEVA) to describe the geometry 
of the vortices that interact with the empennage surface. The equations 
are described in the rotor tip path plane coordinate system. The rotor 
hub location is the origin, the x-axis is positive aft, the y-axis is 
positive to the right (i|; = 90 degrees) and z-axis is pointed up (right- 
handed cartesian coordinate system). The azimuth angle \p is measured 
with respect to positive x-axis, rotating counterclockwise. 

The coordinates of any point n on the vortex trailed from the blade 
radial location r^ (nondimensionalized to R) can be described as follows: 


with 


where 


Alp 


n 




imom 


\t ° ^TPP ((|>n ) (A-1) 


^nt = 


sin (4>n) 

(A-2) 

\t = 

X Y 

TPP a 

(A-3) 

<^n = 


- (n-1) * Ai|; 

(A-4) 

= 

a 

^b 

- ^n 

(A-5) 

^Tpp“ 

Vo 

cos (a^pp)/fiR 

(A-6) 

Xjpp = 

= - 

Vimom^^^ (oi^pp)/fiR 

(A-7) 


is the incremental ip to describe various segments, 
varies from 1 to the maximum number of points chosen, 
describes the current position of reference blade, 
represents age of the wake point, and 
is the uniform momentum induced velocity 
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APPENDIX B 


The mathematical model for the determination of the viscous core radius 
of the blade tip vortices is described in this appendix. In the RIEVA analysis 
the computation of the core radii is required only when the interacting vortices 
come very close to the empennage surface. It is assumed that after roll up, 
the viscous cores can be represented by a constant turbulent eddy viscosity 
model. Under this assumption the core radius r^,, increases with time as given 
by the following relationship: 

r^ = r^ + 3^Yt ('•'a"'*'ru) 
c CO a 

where 

^co core radius at time of roll up, 4^^^, 

Yj. = turbulent eddy viscosity of core, 

3^ = a constant, ( 0 . 716 ) . 

)|/^ = age of the vortex element, 

n = rotational speed 

The magnitudes of r^^, 4'ru \ obtained empirically and input 

to the RIEVA program. 
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EDGE 

VORTEX 

Figure 1. Wing/Vorfex Geometry and Comparison of Measured and Predicted 
Chordwise Pressure Distribution from Reference 3 



Figure 2. Schematic of Main Rotor Tip Vortex/Empennage interaction 
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BOUND 



Figure 4. Modeling of Wing and Its Wake 
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Figure 5. Simplified Model Showing Interaction Effects 
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Figure 7. Flow Chart for Rotor/Tail Vibratory Excitation 
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Figure 8. Vortex Geometry Modification 




Figure 9. Maximum Fb Variation Around Azimuth 
(feet square/sec) 
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BOUNDARY BETWEEN THE TWO PARTS 



Figure 10. An Example of Splitting of Rotor Wake Into Two Parts 
M=0.4, Only Tip Vortices Shown for i/'=30 deg 





TIME (S2t) IN deg 


Figure 11. Variation of Induced Veiocity at Stabiizier 
(At Quarter Chord Point on Wing Centerline) 









Figure 13. Vortex Geometry at Step No. 7 
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Figure 14. Vortex Geometry at Step No. 12 









Figure 15. Stabilizer Airloads Time History 
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